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In this letter, we address the relationship between the statistical fluctuations of grain displace¬ 
ments for a full quasistatic plane shear experiment, and the corresponding anomalous diffusion 
exponent, a. We experimentally validate a particular case of the so-called Tsallis-Bukman scal¬ 
ing law, a = 2/(3 — q), where q is obtained by fitting the probability density function (PDF) of 
the measured fluctuations with a g-Gaussian distribution, and the diffusion exponent is measured 
independently during the experiment. Applying an original technique, we are able to evince a tran¬ 
sition from an anomalous diffusion regime to a Brownian behavior as a function of the length of 
the strain-window used to calculate the displacements of grains in experiments. The outstanding 
conformity of fitting curves to a massive amount of experimental data shows a clear broadening of 
the fluctuation PDFs as the length of the strain-window decreases, and an increment in the value of 
the diffusion exponent - anomalous diffusion. Regardless of the size of the strain-window considered 
in the measurements, we show that the Tsallis-Bukman scaling law remains valid, which is the first 
experimental verification of this relationship for a classical system at different diffusion regimes. 

We also note that the spatial correlations show marked similarities to the turbulence in fluids, a 
promising indication that this type of analysis can be used to explore the origins of the macroscopic 
friction in confined granular materials. 

PACS numbers: 83.80.Fg, 45.70.Mg, 81.05.Rm 


Turbulence is one of the most complex, but ubiquitous, 
phenomena observed in Nature and it is related with the 
underlying mechanisms responsible for the micro-macro 
upscale causing wide-ranging effects on classical systems, 
like macroscopic friction in granular solids or turbulent 
flow regime in fluids [l|-Q. The presence of multiple 
scales in time and space is an additional defy to a com¬ 
prehensive theoretical description, and a particular ef¬ 
fort is made in the literature to perform experiments and 
simulations in order to validate the proposed theoreti¬ 
cal descriptions, particularly Tsallis nonextensive (NE) 
statistical mechanics M- 

A paradigmatic work relating anomalous diffusion and 
turbulent-like behavior in confined granular media was 
presented by Radjai and Roux [l[ , using numerical sim¬ 
ulations, and confirmed qualitatively by experiments by 
Combe and collaborators [3, [s]. Radjai and Roux coined 
a new expression to characterize the analogies between 
fluctuations of particle velocities in quasistatic granular 
flows and the velocity fields observed in turbulent fluid 
flow in high Reynolds number regime, the “granulence”. 
Most of the evidences of the granulence are based in sim¬ 
ulations using discrete element method (DEM) but, un¬ 
fortunately, one can verify a lack of quantitative experi¬ 
mental verification in the last years, limiting the knowl¬ 
edge of the micromechanics of this system based almost 
exclusively on numerical evidences. 


In the present work, we aim exactly to fill this gap 
by contributing with the experimental validation of the 
results obtained by DEM. Specifically, we seek to exam¬ 
ine the findings revealed by Radjai and Roux |kj in a 
detailed fashion, extending the previous works [3, [§] to 
explore quantitatively the relationship between the PDF 
of the velocity fluctuations and the diffusion features of 
the grains. We follow a detailed theoretical description 
for the anomalous diffusion in the presence of external 
driving 1,0- Particularly, a relation between the q- 
Gaussian value from the PDF of fluctuations and the 
diffusion exponent was proposed, which is validated ex¬ 
perimentally here for the first time for a large range of the 
control parameter, differently of previous works where 
this relation was tested only for a single point [ljJ, [l2|. 

In this work, we aim to advance in the route opened 
by Radjai and Roux [flj with three basic goals: 

(z) Explore the low inertial number limit. The inertial 
number I |13l ] measures the ratio between inertial and 
confining forces, from the quasistatic regime (small 
values) to the dynamic regime (large ones) 0 - We 
would like to check if the granulence features are still 
observed in a better stablished quasistatic situation, z.e. 
the experimental one which involves inertial numbers 
around four orders of magnitude smaller than the 
currently reported in simulations 0. 

(zz) Point out the origins of the macroscopic friction. 
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We take advantage of the really quasistatic feature 
of the experimental data to explore the origins of the 
underlying mechanisms of granulence. Here, unlike 
fluid flow, the rigid particles can not fly freely since the 
motion of each particle is hampered by the presence of 
the other particles, and depends on the motion of its 
neighbors. This makes the straining in part controlled by 
geometric exclusions at the particle scale, preventing the 
development of a uniform straining in a sustainable way. 
As shown in (7f, at the limit of large strain-windows, 
it is possible to observe turbulent-like vortexes in the 
fluctuation field which turn out to be associated with 
the energy dissipation and macroscopic friction 16j, 17|. 
(in) Evince the nonextensive nature of the displace¬ 
ment fluctuations. In order to quantitatively analyze 
the data, we have used the Tsallis NE statistical 
mechanics approach. In this context, the PDF of 
displacement fluctuations is not expected to follow the 
normal Gaussian distribution, as in the case of classical 
Maxwell-Boltzmann distribution in thermodynamics. 
In granular systems under loading, the force chains 
engaged along the entire system are a clear evidence of 
long-range interactions [18|. These chains connect the 
microscopic contact forces with the global resistance 
to external forces, as shear for example 0 - Thus, 
it is natural to associate the emergence of these force 
chains at mesoscopic scales with the departure from 
the classical Boltzmann-Gibbs (BG) statistics in these 
systems. 


In our experiment, we have foreseen the possibility 
to quantify the degree of nonextensivity using, the q- 
Gaussian fit of the PDF obtained experimentally [7| . The 
striking accordance observed on the fitted curves, and 
the dependence observed of q as a function of the strain- 
window used to calculate the fluctuations, according the 
reasoning presented here, corroborates the application of 
the NE statistical mechanics on these systems, opening 
an alternative approach to treat these systems quantita¬ 
tively. Besides, by measuring the diffusion of the particles 
along the complete shear test, at different strain-window, 
we are able to associate the g-value measured from the 
fluctuation PDFs with the diffusion exponent a. This is 
a particular case of the Tsallis-Bukman scaling law @, 

( 1 ) 

3 -q 

which can be obtained from the so-called porous media 
equation 0 , a generalization of the classical diffusion 
equation where the linear dependence between the vari¬ 
ance and time is no longer observed 0]: 


dp(x, t) _ d 2 \p(x, t)} 2 9 
dt 9 8x 2 


For a Dirac delta initial condition, the solution reads as 


P q (x , t) = 


1 — 2 — 
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where e q (x) is called ^-exponential, and A q is a constant 
which depends on q and Gamma-function B 0 - 

Equation [3] is known as the g-Gaussian distribution, 
and was used to fit the PDF of displacement fluctuations 
obtained experimentally. Figure [2] shows the results for 
the PDF of fluctuations and the corresponding fit func¬ 
tion at two extremal values of Ay considered in the image 
analysis. 

In the case of anomalous diffusion, it is shown that the 
variance follows a power law with time: 


(x 2 ) oc t a = (x 2 ) <x t 3 ~i , (4) 


where a is the diffusion exponent equivalently expressed 
as a function of q by using Eq. [TJ Note that, in Eqs[2]to 
Q] x stands for a fluctuation of displacement as we will 
see below. It is interesting to observe two special cases: 
when q = 1 , the variance is proportional to time which 
corresponds to the normal diffusion behavior; when q = 
2 , the ballistic diffusion limit is reached. At intermediate 
values, we get large distributions with marked tails. The 
variance, calculated as the time-integral of p q , diverges 
for q > 5/3, and converges otherwise. Thus, if several 
independent convolutions are applied, p q approaches a 
Gaussian distribution if q < 5/3, and it approaches a 
Levy distribution for q > 5/3 [ 6 j. 

We performed quasistatic simple shear tests with ly2e 
apparatus which is fully described in (21, 0. The gran¬ 
ular packing is made of pilings of cylindrical rods which 
mimics a 2D granular material enclosed by a rectangular 
frame, with initial dimensions of 0.56 m x 0.47 m. Then, 
the vertical sides of this rectangular parallelogram are 
shortened or elongated to apply a constant normal stress 
in the vertical direction, a n = 50 kPa. These two vertical 
sides are tilted up to 7 = 15° while the two other sides are 
kept horizontal with a constant length - Fig. [I] The pack¬ 
ing was made of 5471 wooden rollers (6 cm long) with ten 
different diameters ranging from 3 mm to 30 mm, ap¬ 
proaching to a uniform distribution. To ensure a quasi¬ 
static transformation of the sample, it is sheared very 
slowly - the corresponding shear rate 7 is 4.5 x 10 ~ 5 s" 1 . 
This ensure a very small inertial number 0 (J = 10 " 9 ) 
when compared to what is applied in DEM simulations 
(10 -3 to 10 -5 in the best cases) ,15, 23]. During the test, 
kinematics of grains are measured by means of Digital 
Image Correlation (DIC) jH,0 from 80 MPixels digital 
images of the sample where rollers look like disks, Fig. |T] 
A specific DIC computer program was developed to track 
rollers here assumed as rigid bodies @ which allowed a 
sub-pixel kinematics measurement to track grains with 
an error of ±0.05 pixels [26 ]. 
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FIG. 1. The fluctuating part of the rod displacement u(y, A 7 ) 
over-plotted on the corresponding digital image, obtained 
from DIC technique. 7 = A 7 = 0.1. Inset: a detailed view 
of the speckled rods. 7 and cr n are the shear angle and the 
vertical stress imposed, respectively. The shear strength is 
measured all along the shear test. 


At the macroscopic level (sample scale), the stress- 
strain curve measured during the shear test exhibited 
hardening up to 7 « 0.06 followed by softening until the 
end of the test (curve shown in [§| as well as several other 
mechanical properties like peak stress ratio, macroscopic 
friction angle etc.). Pictures were shoot every St = 5 s 
throughout the test, corresponding to a shear strain in¬ 
crement A 7 = St 7 ~ 2.4 x 10 ” 4 between each shot. To 
assess the displacement fluctuations, we consider two dis¬ 
placements of each particle during a shear increment A 7 . 
The first is the actual displacement Sr( 7 , Ay) from 7 to 
7 + Ay. The second displacement, Sr*( 7 , Ay), is ficti¬ 
tious and corresponds to an affine motion resulting from 
an homogeneous straining at 7 and during the shear- 
increment Ay. It is assessed from the motion of the four 
rigid sides of the apparatus ly2e. With these definitions, 
the fluctuating part of the displacement is the difference 
between the actual and affine displacements. Thus, the 
normalized displacement fluctuation u(y, Ay) is defined 
by: 


u(y, Ay) 


[<5r(y, Ay) — <5r*(y, Ay)] /d 

Ay 


(5) 


where d is the mean diameter of the rollers. One may no¬ 
tice that the normalized fluctuations can be interpreted 
as a local strain (grain scale - numerator) compared to 
the global strain (sample scale - denominator). 

A displacement fluctuation field v is plotted in Fig. |T| 
for a given shear increment. We notice an organization in 
structures like vortexes that reminds the ones of turbu¬ 
lence in fluids. These structures found their origins in the 



FIG. 2. Probability density functions of the horizontal com¬ 
ponents of the fluctuating displacements tracked during two 
different increments of shear strain (Ay = 7.3 x 10 ~ 4 and 
Ay = IIP 1 . The scatters correspond to experimental data, 
and the solid lines correspond to regressions of function p q 
(Eq. K3) that allow for the assessment the q-values. 


rearrangement mechanism of the grains, since that the 
elements interfere with each other in their affine move¬ 
ment. This is, in other words, the deviation from the 
affine field due to steric exclusion forming pat terns ob¬ 
served with discrete element modeling [13j, |27| [28| and 
more rarely in experiments [17] ■ Their dynamics depend 
both on y and Ay, evolving gently under shear when 
Ay is large (> 0.04) and very rapidly for small values 
(Ay ~ 2.4 • 10~ 4 ). The characteristic lengths depend 
strongly on Ay, with vortexes of a few tenths of grain 
mean diameter for large values of Ay, and, to the con¬ 
trary, for small values of Ay these structures are not well 
defined, and long range correlations are observed @. 

The PDFs of the horizontal component magnitude of 
normalized displacement fluctuations are shown in Fig. [2] 
for two different increments of shear strain: Ay = 7.3 • 
1CP 4 and Ay = 10 _1 . We observe a broadening of the 
PDF from a nearly Gaussian distribution (for large Ay) 
to a wider distribution (for small Ay). 

The dependence of the q —exponent with the strain- 
window used to calculate the fluctuations is shown in 
Fig- El for experimental and simulation data. Two re¬ 
markable features can be observed in this plot: first, in 
the limit of large strain-window, when the abscissa goes 
to zero, q —> 1 , indicating the limit when normal diffu¬ 
sion and the BG statistics are satisfied. Note that it is 
possible to test larger values of Ay in DEM simulations 
which confirm the limit q —> 1 (data shown in the inset 
of Fig. [3]). This is exactly what we expect for this limit, 
once that the particles typically experience several colli¬ 
sions and rearrangements, approaching to the molecular 
chaos hypothesis. 

In the other limit, for vanishing strain-window, the q- 
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FIG. 3. Evolution of the measured g-value as a function of 
squared inverse of the strain increment for the experiments 
and simulations. The dashed line corresponds to a regression 
with the function g(l/Ay) = 1 + atanh(b/ Ay), with a = 
0.521, b = 0.096. Inset: the same plot for data from a DEM 
simulation that highlights the limit q = 1 when Ay —¥ oo. The 
fitted parameters for simulations were a = 0.387, b = 0.057. 


value attains a plateau, with q ~ 3/2. This observation 
can be interpreted as a sign of the long range correlations 
imposed by the force chains at this short time scale. Once 
the value measured for q in this limit is lower than 5/3, 
one can expect that for large strain-windows a Gaussian 
distribution would be recovered, since it correspond to 
successive independent convolutions of g-Gaussian dis¬ 
tributions. 

These features were observed both in experiments and 
simulations, no matter the differences among the systems 
(periodic boundaries in horizontal direction in simula¬ 
tions, different number of particles and inertial numbers 
etc), proving the robustness of the result. 

Analyzing the results as a whole, we can sketch a phe¬ 
nomenological scenario to explain the observations: in 
the limit of large Ay, we observe a tendency to agree with 
the BG statistics, with q —> 1. This limit corresponding 
to the transition from meso- to macroscopic scales, and 
we observe the formation of vortexes in the spatial dis¬ 
tribution of fluctuations, as evinced by Fig. [Q These 
vortexes, with few grains diameters in size, interact each 
other to dissipate the excess of energy due to external 
loading, in analogy with the role of the vortexes in tur¬ 
bulent flow 0. The nature of the interactions of these 
structures is purely stochastic, which acts as a precursor 
for the macroscopic friction. The broadening of the dis¬ 
placements fluctuations distribution is usually attributed 
to the energy cascade from larger to lower scales, that is, 
from large vortexes to the small ones 0. 

On the other hand, we have g ~ 3/2 for vanishing 
strain-windows limit, Ay —>- 0. This result indicates the 
presence of long-range interactions and anomalous diffu¬ 
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FIG. 4. Verification of the Tsallis-Bukman scaling law for 
different regimes of diffusion, (top) Evolution of the measured 
diffusion exponent a as a function of 1/Ay the dashed line 
is a direct application of the scaling law from the fit of the 
values shown in Fig. 02 ap(l/Ay) = 2/[3 — g(l/Ay)], where P 
stands for predicted. (Inset) a selected curve diffusion where 
x is the displacement fluctuations; it allows the assessment of 
the diffusion exponent a. (Bottom) Measure of the deviation 
of the data relative to the scaling law prediction, as a function 
of 1/Ay, showing a agreement on the order of ±2%. 


sion. Considering the absence of spatial structures on the 
fluctuation field, it is clear that this limit is dominated 
by the force chain dynamics. Force chains can span all 
along the system, but are very fragile, implying short life 
times. The displacement of grains belonging to a force 
chain are strong correlated spatially, but this correlation 
is not verified for temporal scales. 

Thus, we can conclude that the window used to mea¬ 
sure the PDF particle displacement fluctuations in the 
system plays a crucial role in the statistics that will be 
obtained. Basically, it is possible to explore the micro¬ 
macro transition on the PDF distributions, from a corre¬ 
lated regime dominated by the force chains to a frictional 
stochastic one, dominated by spatial vortex interactions. 
This conclusion has a striking implication for any anal¬ 
ysis concerning the measuring of displacement fluctua¬ 
tions, since it unveils how the observation procedure can 
alter the conclusions even in a relatively simple diffusion 
experiment. 

To quantify the diffusion of the grains along the com¬ 
plete shear test we basically computed the average dis¬ 
placement of each grain as a function of time (shear in¬ 
crement y), but with different sampling frequencies de- 
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termined by the strain-window A 7 . Following the rea¬ 
soning presented above, and Eq.[4j we should expect two 
extreme regimes for the diffusion, analogously to the ob¬ 
served for the g-value: an anomalous diffusion regime 
with a ~ 4/3 for vanishing strain-window, and an asymp¬ 
totic regime with a —> 1 for large shear increments. This 
is indeed what we can observe in Fig. [4j where we have 
verified the Tsallis-Bukman scaling law (Eq. [U. It is 
important to stress that the dashed line in Fig. [4] is not 
a direct fit, but rather the curve obtained in Fig. fusing 
the Tsallis-Bukman scaling law. To our knowledge, is the 
first time that this relation is verified for different regimes 
of diffusion. This striking result reinforces the use of the 
Tsallis NE statistical mechanics to describe strong corre¬ 
lated systems, as in the case of confined granular material 
under shearing. 
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